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ABSTRACT 


A discrete potential element approach to subsonic numer- 
ical lifting surface theory has been developed and shown to 
be practical in predicting the nonsteady loading on 
harmonically oscillating, medium to low aspect ratio wings. 
A unique method of including the wake effect in the wing 
Remnel function matrix prior to solution of the singular 
integral downwash equation was devised, thus greatly 
simplifying the velocity potential formulation. In addi- 
tion, termination of the effective wake a finite distance 
downstream of the wing was investigated, with wing loading 
Pound FO converge Co within one per cent in an effective 
Weieow length of four root chords. 

oomorsenelomclcmMent method has also been extended to 
the case of an oscillating wing, cantilevered from a cylin- 
drical fuselage, to investigate nonplanar interference 
effects. This interference in wing loading, while of rela- 
tively small magnitude, does exist in both pressure amplitude 
and phase angle distributions, and is, therefore, of impor- 
tance in three-dimensional stability analyses of wing/body 


CO 1Suravions . 
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ieee RODUCTION 


iiemecrodynamie analysis of the forces over an oscillat- 
ing wing, situated in a steady flow, is the basic problem 
Censidered by nonsteady lifting surface theory. The deter- 
mination of these nonsteady pressure distributions is of 
perme importance in the consideration of stability. A basic 
objective of flutter analysis, for example, is to determine 
the aerodynamic loading on a wing oscillating with small 
amplitudes in a definite mode and.with a given frequency. 
Thus the problem is formulated on the basis that the motion 
and deformation of the lifting surface are either known or 
composed of known elemental modes, and is amenable to har- 
monic modal analysis. A downwash integral equation, based 
on the transformed wave equation for linearized potential 
flow, relates loading (either in terms of velocity potential 
or pressure) over the oscillating wing to the normal velocity 
induced by the wing motion, and must be solved numerically 
Meomoobcain this wing loading. 

The analysis of this paper will deal primarily with the 
subsonic case, this being the one in which the largest 
number of applications of lifting surface theory to engineer- 
ing problems are to be found [1]. Most of the work in sub- 
sonic nonsteady lifting surface theory is based on the 
development of the general kernel function approach by 


Poaemeieim 1940 [2], which was first formulated successfully 





for the computer by Watkins, Runyan, and Woolston in 1955 
[3]. More recent techniques have been developed by Stark in 
1964 C4] eMemeacehitauim 1963 15]. Solution of the singular 
downwash integral equation by the usual lifting surface 
methods is accomplished through assuming the loading to be 

a series of preselected functions with unknown coefficients, 
and then determining these by satisfying the normal velocity 
eemalvion with some sort of collocation technique. The 
downwash and the loading are related by kernel functions, 
which are themselves singular, and which must be numerically 
evaluated. 

Mimapitye Of the fact that lifting surface theory has 
enjoyed a great deal of success, certain shortcomings, 
discussed in Refs. 1, 6, and 7, do exist: 

(1) Lifting surface theory is not as computationally 
economic as steady or two-dimensional methods. 

(ii) Wing loading functions, which incorporate the 
proper singular behavior, have to be assumed a priori. 

(144i) Uncertainty exists as to the method of locating 
Somlocation points. 

(iv) Techniques for analyzing control surface effects 
have not been adequately developed, partly because of (ii). 

(Vee Lifting surface theory has been primarily restricted 


woeplenar surfaces. 


In order to attack these problems, Houbolt [6] proposed that 
the wing loading be represented by concentrated pressure 
loads in the form of Dirac delta functions. Using this 


loading representation, the kernel function and downwash 


1S: 





equation were formulated and analyzed for certain restricted 
Gases. However, no general numerical solution of the lifting 
surface problem was developed. 

The present study is concerned with developing such a 
method of solution to the subsonic, nonsteady lifting surface 
problem, in terms of discrete loading elements. However in 
this case wing loading is expressed in terms of velocity 
potential, rather than pressure, since this formulation 
allows a more direct approach to nonplanar configurations. 
Havitand [8] reports using a similar approach to the steady 
planar case with consistent results. 

iicmaiimOonwvnis discrete potential element method is to 
provide a new approach to nonsteady, subsonic lifting 
surface theory which eliminates, or minimizes, the problems 
previously indicated. Theoretical development of this 
aem@eacn was made from basic principles, and the resulting 
formulation applied in a computer program which analyzes 
Memmemlcally Oscillating wings in subsonic flow. A unique 
way of including the wake effect was also developed which 
fieauly simplified the velocity potential solution of the 
numerical problem. As an extension of the theory, and an 
indication of the versatility of this discrete element 
approach, a second computer program was developed to analyze 
interference effects on oscillating wings mid-mounted from 
steady bodies. To the author's knowledge, this type of 
wing/body interference investigation has not previously been 


presented for the nonsteady case. 


eG a 


eS ie 





it. GENERAL PROBLEM DEVELOPMENT 


Following the method of development by Garrick [9], the 
governing equations and boundary conditions are linearized 


Byeapplication of small perturbation theory. 


A. POTENTIAL FLOW EQUATIONS 
The governing equations for adiabatic, irrotational, 


inviscid flow, as developed in Ref. 10 are: 


omc rnuit y 
oiey otl= 0 (2 aay) 
ot 
Momentum , 
30 , 2 1 
apt VU = = | VP (2a) 
Energy 
= = constant (y the specific heat ratio) (20) 
p 


Since the flow is irrotational, the curl of the velocity 
vector vanishes (VxU=0) and a velocity potential function 


Cam be defined, such that 


=e Gon 


cy 


If further, a perfect gas is assumed, then from the energy 
equation the well known relationship 


dP/dp = ac C225) 


is obtained for constant entropy flows. 
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The continuity and momentum equations can now be rewritten as 


i Doe -6 = 
Cebe + Vo = 0 (2.6) 
2 2 
d® 8 as aa a_ 
v( 38 aa fara: C2 
De, : : ; ; 3 es 
where Dt iS the substantial derivative defined by [asi vy. 


Combining equations (2.6) and (2.7), a single equation 


in @ is obtained. 


ve - 3 (24 t-9]% = 0 


ac ot 
Or 
o 
veo - = PS = 0 (2.8) 
a Dt 


Sess 2S the governing partial differential equation for 
general nonsteady, nonviscous, potential flow of a perfect 
Peaomeeequation (2.6) is not limited to small disturbances, 
and has the form of the classic wave equation in terms of 
Moma ubstantial, or convective derivative. Thus, the flow 
disturbance represented by the velocity potential is con- 
yected by the local fluid, or stream velocity, and propogated 
as a wave which spreads at a rate equal to the local speed 
Seo ound 

In this general form, equation (2.8) is highly nonlinear 
Siemoecem re che quadratic terms in the convective operator 
Uv, and to the interdependence of the velocity potential 
Ga@menae local speed of sound. Boundary conditions, in 


general, will depend on the location and form of the moving 


18 = 





Peay, O11 Ciscontinuities along streamlines, on shock waves, 
and on flow conditions at infinity. There exist no general 
methods for obtaining solutions to equation (2.8), and it 
is necessary to consider small perturbations, to linearize, 
moereduce the number of equations, or to try schemes of 


successive approximations to attack the general problem. 


B. OMALL PERTURBATION THEORY 

Nesey Of the theoretical developments in aerodynamics, 
including nonsteady aerodynamics, are based on the concept 
Of small perturbations. There are two aspects to the prob- 
Hemeot tlow past a body creating small disturbances in an 
undisturbed mainstream, linearization of the governing 
differential equations and further linearization of the 
Beundary conditions. 

Manearization of the governing equations requires the 
problem to be posed with a compressible fluid in a uniform 
stream flowing with velocity U, in the positive x direction. 


Smalt disturbances in the flow are defined by 


U = (Uta) 4+ V4 + wk = Uli t+ 

p=p.t+op d= >, + $ 
= => _ 

eect: 2 with u = Vo 


The perturbation velocities are considered small with respect 
weetesse., and U_-a_, thus ruling CuGeoreonsomle 1 low. ihe 
perturbation density and pressure are also restricted, so 


chat 
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p/p, and p/p, «l 


femoidering © as the perturbation velocity potential, 
the governing equation (2.8) becomes 


Z 


Oat 0 (25S) 


O= ab re) Pe) 
%3-Sik+ud 





where higher order terms in the UV operator have been 
neglected. This linearized governing equation may now be 
applied to flow problems within the initial small perturba- 
tion assumptions. Although equation (2.9) was developed for 
a uniform flow in the positive x direction relative to a 
meee —1ixed coordinate system, it also applies to the case 
fea wing moving with uniform velocity U, in the negative x 
Garection with the coordinate system fixed to the wing. 
Within the concept of small perturbations, the momentum 


equation becomes 


> 
ohl a 1 o- 
+uU, = - => Ww (2710) 


BUesebivuving Che perturbation velocity potential (i =VS) 
into this equation, a linear relationship between potential 


‘and pressure is obtained. 


2 Eee 


Therefore p satisfies the same differential equation (2.9) 


as does 6. 


If harmonic motion were considered, the perturbation 


variables would have the form 


20 





iwt - ae 


Wear eee uets) 0 P(X,y,z)e 


eo 
it 


M@eyee eo”. etc. 


ct 
Il 


and the governing equation would become 


Cie 6 = 0 (2.12) 


with the relationship between pressure and velocity potential 
is 3 
soe Warn 7 10) o (2s) 


In the above equations the time dependence aoe has been 
factored out. This relationship between p and @$ can be 


integrated to give 


oe = iwé/U, 
Ors sz ) =e ape f Cary. ae ae (2514) 


co 4 


where far ahead of the disturbance, the perturbation poten- 


Tial is assumed to be Zero. 


C. LINEARIZED BOUNDARY CONDITIONS 

The small disturbance assumption, resulting in the 
linearized governing differential equation (2.9), implies 
aieleadevyiations from the uniform flow. Therefore it is 
also necessary to linearize the boundary conditions for the 
physical problems considered, such as properly oriented thin 
wings and slender bodies. These boundary conditions consist 
of the mathematical formulation of several physical and phe- 
nomenological statements which are normally grouped as 


moOllows: 


eal 





(1) Surface boundary conditions. No flow can pass 
memouehn the wing or body surface, that is to say the total 
Mee 1S tangential to the body surface. 

(ii) Edge conditions. Sufficient viscosity is present 
in the "nonviscous" flow to determine the flow pattern near 
sharp edges. Thus in subsonic flow the well known Kutta 
Condition can be applied to the wing trailing edge, requiring 
that the surface pressure difference remains continuous 
near, and vanishes at, the trailing edge. A similar condi- 
tion should be satisfied at the side edges. Through the 
method of matched asymptotic expansions, Landahl [11] veri- 
fmeed that at trailing and side edges the pressure difference 
veer approach zero with infinite slope due to weak singu- 
larities in the first derivative. Conditions at the leading 
Saremnave not been fully explored, but for small disturbances 
memrounded leading edges it is sufficient to require that 
the total integrated force be finite and the singularities 
in the pressure distribution be of the proper order [9]. 
woes latter requirement is fulfilled with the pressure 
varying as the inverse square root of the distance downstream 
of the leading edge [12]. 

(iii) Wake conditions. The free vorticity shed from 
Pace trailing edge in the wake is such that its circulation 
Boeether with the bound circulation vanishes in accordance 
with the Helmholtz-Kelvin Theorem. It is assumed that the 
shed wake remains where it is formed, floats without mutual 


interference along streamlines, and forms a continuous sheet 


ee 





memearocOncinuity coplanar with the nasleaye iaewimiendirection of 
Mmmene. Hdge effects and roll-up of the sheet at infinity 
are ignored. 

(iv) Conditions at infinity. The flow is considered 
uniform at infinity, and perturbation waves are required to 
Mroepagzace away from the sources of disturbances and to 
Memave properly at infinity. 

ijeaddition FO the above, for supersonic flow account 
must be taken of zones of influence and dependence, as in 
the method of characteristics. In supersonic flow there are 
eoeorouner problems which do not arise in the subsonic 
analysis, such as the difference between subsonic and super- 
sonic edges. 

To formulate the main boundary conditions, a thin wing 
as considered lying in the x-y plane, creating small dis- 
Pemeaces in a uniform stream U.. flowing in the positive x 
direction. The linearized tangential flow condition on the 
wing surface Be takes the form 


OZ OZ 


w(X,y,Z.,t) = a af ie Tae C25 15s) 


The normal fluid velocity w can be expanded in a power series 


eeeue Z = 0, so that 


= = dw = 
w(X,y,Z.,t) = le we AIG tz)! ae 4 os 
Z=0 
Sipe 
i He fed ze aaa (Za) 
' Moz iz=0 





mmc bon, problems of camber and thickness can be separ- 
ated from wing motion and angle of attack in the linearized 
Memmulation. Therefore, neglecting higher order terms in 


equation (2.16), the tangential flow condition becomes 


7 Z az, 
w(x,y,0,t) = —~ + U aoa (2347) 


co 


Q}; @Q 


wee z(x,y ,0,t) represents the mean camber surface posi- 
een Of the wing. 

Since w is an even function of z (its value is unchanged 
for the top or the bottom side of the mean surface), the 
velocity potential b is an odd function of Zz, as is the 
Perturbation pressure p through the linear relationship of 
equation (2.11). Therefore, in the plane of the wing b 
equals zero outside of the wing and wake. 4 does not equal 
Zero in the wake because of the discontinuity in the u 
component of the perturbation velocity across the wake. The 
perturbation pressure is also zero in the x-y plane, except 
at the lifting surface where the pressure difference Ap 


between the top and bottom surfaces is given by 


DieGsy 50-55) re (Oe a Chey B, 


AND Coe aie ce) 


—AgiGene0.t) = —20 ols t sak C26) 


Figure 1 summarizes the formulation of these boundary condi- 
tions for a lifting surface in the x-y plane. 

These linearized boundary conditions along with the 
linearized governing differential equation (2.9) can be used 


to attack a wide range of perturbation problems in subsonic 
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FIGURE L 
COORDINATE SYSTEM AND BOUNDARY CONDITIONS 
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Ze 





and supersonic flow. Miles [13] has shown that a sufficient 
Bemoiti10On for linearization of the equations for the pertur- 


Beeton potential is any one or combination of the factors 


IM? - 1] > 0 (6°73) 


k > 0 (62/3) 


il W3 
AR? 9 (6 ) 


ere oO, M 6, ké, and kM_5 must be sufficiently small. 
The perturbation equation is nonlinear only when the 


femebowing conditions are jointly satisfied 


hee tf = 0 (6°/3) 


1 _ 1A ° 
a = 0 (6 ) 


6 is here intended to be a thickness, angle of attack, camber, 
or motion parameter serving as a measure of the flow 


disturbance. 
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ieee PNG SURFACE THEORY 


fee BASIC APPROACH 

The basic nonsteady problem considered in lifting surface 
theory is that of a linearized thin wing harmonically 
Zeemelacing in the normal direction in some prescribed 
deflection mode, creating small disturbances in the uniform 
mainstream. Particular solutions to the linearized governing 
differential equation (2.9) are required which give pertur- 
bation velocity potential or perturbation pressure distri- 
butions on the wing, and which satisfy the essential boundary 
Pemearcions discussed in Section II1-C. Such a direct 
Semuclon has been possible for only a few special cases, 
such as two-dimensional incompressible flow as in Refs. 14 
and 15. For the more general three-dimensional case, an 
integral equation formulation must be employed which makes 
use of a downwash kernel function, analogous to an influence 
coefficient development. 

Considering the time dependence ge Pave shia mond c 
analysis as being factored out, the integral equation relating 


the downwash amplitude to the wing loading has the form 
w(x,y,0) = ae ieee 0) KC xe 5 ym 50 )dédn C310 

The loading L can either be the velocity potential or the 

pressure distribution over the surface, and K is the kernel 


function which denotes the normal velocity (downwash) on the 


Pamemacex,y,0) due tO a Unit acoustic singularity located 


oa 





at (&,n,0). In other words, the surface is modeled by a 
meer ouLion Of radiators, formulated in terms of either 
Meeetivbial or pressure which satisfy the linearized governing 


equation 


es 2 n> Ae (2.12) 


w is the downwash at point (x,y,0) on the wing, determined 
DY the harmonic motion of the surface in fulfilling the 
tangential flow requirement. <All of these quantities are 
Gomplex in the spatial domain due to the nonsteady nature of 
the problem. 

Employing the velocity potential approach, the kernel 


miretb10n is defined by 
K = w(x,y,0) = 96/9z : Cac) 


Wigeme > 1S the solution to equation (2.12) for an acoustic 
Mmeretavor Of unit strength. In this formulation, the inte- 
gration (3.1) must be carried out over the surface of both 
wing and wake, since the potential does not go to zero in 
the wake region. This complicates the integration process 
and has proven a major disadvantage of the potential approach. 
If the wing is modeled by a pressure distribution, the 


kernel function is denoted by 


—-iWXx ee 
Sees 2 ee Yoo 
K=w(x,y,0) =———- Lim ~/fS p(&,y,z)e aé 
oo Ue 7>O0 dz =O 


(373) 
from equation (2.14) where p is now the solution to equation 


meee tor an acoustic radiator of unit strength. The 
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mmieeeracloOn here need only be carried out over the wing 
surface, since the pressure goes to Zero in the wake. 
Because of this, methods based on the pressure formulation 
have been more widely employed. This kernel function (3.3) 
is more complicated than that developed from equation (3.2) 
and has not been integrated in closed form due to singulari- 
mess at (y-n)=0 and (x-&)> 0. However these singularities 
have been isolated and expressed in forms which can be 
handled by numerical procedures [3]. 

The acoustic singularities used to model the wing surface 
in subsonic flow are harmonically pulsating doublets (dipoles) 
merenoed in the Z, or wing-normal, direction. Doublets are 
employed in order to represent the pressure difference 
between the two surfaces of the wing, while <n Supers ons 
flow acoustic sources, or monopoles, are used due to the 
independence of the two wing surfaces in that flow regime. 
Kussner [2] developed the general formulation of the pressure 
@eubplet Kernel function, through use of the Lorentz trans- 
MemiiatiOn, in order to apply solutions to the classical wave 


euacton of physics 


to the governing differential equation of lifting surface 
wmeory (2.9). This method of attack is used in Appendix A 


to develop the kernel functions used in this report. 
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B. NUMERICAL SOLUTIONS TO LIFTING SURFACE THEORY 

fo find the potential or pressure loading on the lifting 
surface actually requires the inverse solution of equation 
(3.1), since the downwash velocity distribution is deter- 
mined by the type of surface motion prescribed through 
equation (2.15). The numerical methods that nave been 
developed in subsonic lifting surface theory rest on this 
m@iverse solution to the singular integral equation. The 
usual procedure is to assume the loading to be a series of 
Meeosetrected functions with unknown coefficients and then 
to determine these by satisfying the downwash velocity 
feetribubion exactly in a set of collocation points [16, 17], 
Or approximately in the least squares sense in a larger set 
memcentrol points [4, 18], or by satisfying certain integral 
memacvlons derived from variational procedures [19, 20]. 

Ashley and Landahl [21], and Landahl and Stark [1] have 
presented survey papers summariZing the status of numerical 
difting surface theory, and presenting formulations of the 
integral downwash equation and the kernel function in the 
various methods of attack, such as: 

(i) Velocity potential. As previously discussed, the 
SaegeretUne>ion used in this formulation, while singular, is 
much simpler than that employed with the pressure loading 
approach. However, integration must be carried out over 
both wing and wake. The velocity potential formulation has 
maeeemaloyed by Jones [22] for M, = 0 in 195e. 

(ii) Acceleration potential. Formulating this approach 


in terms of an acceleration potential defined by ee -2y [23], 


30 





allows the use of the same kernel function as in the velocity 
Memeritial formulation, but removes the necessity of integrat- 
ing over the wake, since Cy = Q there. However, the solution 
is not unique, since multiples of eigen solutions with 

ge a O on the wing may be added, and integration must be 
extended into the region ahead of the leading edge to achieve 
uniqueness. This formulation has been used for two- 
emensiOnal cases and also for certain three-dimensional 
cases [24]. 

(iii) Pressure formulation. The original development 
by Kussner [2] has been formulated for the computer by 
Watkins, et al [3] and by Richardson [16]. This procedure 
provides direct determination of the pressure with integra- 
tion required only over the wing surface, but the kernel 
function is highly singular and needs to be evaluated through 
numerical quadrature, thereby increasing computer time 
considerably. 

Various refinements of the basic approaches discussed 
above have been developed, such as the integrated acceler- 
Peron potential by Stark [4] and the advanced velocity 
potential [1], but these do not change the basic problem 
formulation, nor the basic advantages and disadvantages of 
meenwconcinuous loading function methods. 

Dielit cing surface theory the loading is expressed as 
a linear combination of preselected functions, with coeffi- 
cients to be determined by satisfying the tangential flow 


@eneition through the integral equation (4.1). In selecting 
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the loading functions, known results from two-dimensional 
mrecompressible flow theory have been used. The loading is 
mewaracred into chordwise and spanwise functions, so that a 


M@meessure function element would be represented as 


oe f fé)e, (n) 


These loading functions have been represented by Fourier 
expansions [17], power series [25, 26], Tschebycheff 
polynomials [27, 28], Legendre polynomials [4], and so on. 
dhe basic requirement is that these functions themselves 
satisfy the applicable boundary conditions from Section II, 


Meoeas leading, side, and trailing edge behavior. 


fee oULUTTONS BASED ON DISCRETE LOADING LINES 

igesoveady flow, the lifting line theory- of Prandtl 
employed a discrete vortex lifting line, with variable 
strength, placed at the quarter chord line. Wieghardt [29] 
exvended this method for rectangular wings by using several 
meecine® lines, while Rubbert [30], Dulmovits [31], and 
Hedman [32] assumed constant strengths of the lifting lines 
in subintervals in a vortex lattice approach. 

The vortex lifting line in steady flow corresponds in 
the nonsteady case to a potential doublet strip with strength 
varying harmonically in the streamwise direction. The use 
Of a discrete lifting line approach to the nonsteady case 
Heoeiirst proposed by Jones [22] for the incompressible 
case, and by Runyan and Woolston [25] for general subsonic 


iow . 
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Mematetine Steady case, this lifting line approach for 
memovecady motion has been extended to the doublet lattice 
method by Stark [1], and by Albano and Rodden [33]. The 
Gewnwash from the doublet lattice strip is closely related 
tO the kernel function discussed previously in this section. 
pince the surface loading is replaced by discrete doublet 
Ser ipos, it iS not necessary to assume the loading functions, 
@emreguired in the continuous function approach to lifting 
surface theory, and the integral downwash equation (3.1) is 
replaced by a matrix equation relating downwash on the wing 
mo doublet lattice strength. This approach then removes 
one of the shortcomings to general lifting surface theory, 


but still retains the remaining features. 


D. NONPLANAR CONFIGURATIONS ° 

Limited work has been done in the area of general non- 
beater configurations, since interference investigations 
have been primarily involved with combinations of lifting 
Meer esuriaces. The kernel functions for interfering planar 
surfaces, as developed by Davies [34] and by Landahl [35], 
become more complicated but do not contain any new singulari- 
ties. In the selection of suitable loading functions, the 
memevior at the intersection of lifting surfaces must be 
taken into account. The T-tail has been treated by Stark 
[4] and by Davies [35], while calculations for a delta wing 
with folded tips have been given by Vivian and Andrew [36]. 
Lashka [5] has analyzed the effects of wing tip pylons as 


well as interfering planar surfaces [37], while Rodden [38] 


33 





has included wing/body and wing/pylon effects in nonsteady 


wing loading in subsonic flow. 
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hee PrsChe th POTENTIAL ELEMENT DEVELOPMENT 


A. WING ANALYSIS 
ime fSasic Formulation 
The integral equation relating downwash to wing 
loading (either velocity potential or pressure) in lifting 


surface theory, as discussed in Section III, is 


w(x,y,0) = JBE L(&,n,0)K(x-&,y-n, 0) d&dn (3 sail) 
S 


where harmonic motion is assumed so That 


wi 


mex 425 b ) w(x,y,z)e 


aL hie 


L(x,y,z,t) L(x,y,zJe 


<é 


In the subsequent development, the time dependence will be 
memoeaered factored out, so that all symbols will refer to 
complex valued spatial variables. The downwash w is deter- 
Genea by the boundary condition of no flow through the wing 
surface as given for general motion by equation (2.7), and 
Bere for harmonic motion 


OZ 


= fe 7 
W J Gr + iwz Gl a.) 


where Z is the amplitude of the surface motion prescribed 
mereune problem. The form of acoustic radiator used to 
model the wing is taken as a potential doublet, with axis in 
gemzy Or Wing normal, direction. This dipole singularity 


B=wemployed because of the requirement to sustain a pressure 
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aifference across the wing in subsonic flow. K is the kernel 
function relating the downwash at point (x,y,0) on the wing 
Mento a Unit potential doublet located at (&,n,0). 

In normal lifting surface theory, the loading L would be 
represented by continuous chordwise and spanwise functions 
which satisfy the surface boundary conditions developed in 
Section II. However, in the discrete potential element 
approach this loading is represented by a network of Dirac 
delta functions formulated in terms of the perturbation 
Beeeecoy pOLencCial. Thus the velocity potential distribution 
Over the wing is replaced by a series of point functions in 
the form of harmonically oscillating potential doublets with 
ewes in Che Z direction. In sake aoI the wake must also 
contain a network of these doublets since the velocity 
potential is not zero there (Section II) — the integration 
region of equation (3.1) must include both wing and wake. A 
mewoeesentation of this potential Cenpletecria 1S Shown in 
Figure 2. | 

ieecie Giscrete formulation, the integral equation (3.1) 


is replaced by the matrix equation 


wade SUT ce) 


where Ws is the downwash at the yon point on the wing, and 


>, is the strength of the ye 


and wake in the grid of Figure 2. The location of the down- 


potential doublet on the wing 


wash control points are coincident with the potential doublet 


positions at the center of each wing control box. The choice 


Seechis location will be discussed in Section IV.A.4. In- 
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modeling the wing, the velocity potential and the pressure 
are considered constant over each of the control boxes, with 
meeroximations to continuous distributions given by these 
Values at the control box midpoints. 

The kernel function matrix elements represent the down- 
Wash due to potential doublets of unit strength which satisfy 


maemwecoverning linearized differential equation 


ieee = 0 | (2.12) 


As developed in Appendix A, this kernel function has the 


form 


Gd) 











=e ; : Ww) 
a (l+i R)exp{i [M(x,-x,)-R]} (4.3) 
ig haR? a 3° a Be oa 


Lee) co 


where 





2 z 
(x5-X,) + 8 (Fars? 


iomeoenls formulation the kernel function is singular only 
Miemeecne GOwnwash control point is coincident with the poten- 
tial doublet location (i=j). The handling of this singular- 
ity is vital to this development and is discussed in detail 
in Section IV.A.4. 

iafemsalccreve lodding element approach to the nonsteady 
wing problem requires the solution of the system of linear 
complex equations (4.2) to determine the velocity potential 
vector {¢}. The downwash vector {wi} is specified by the 
Wing motion through ane Cris so that the solution has 


the form 
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ee (KI! Cut (4.4) 


where {¢} and {w} cover both wing and wake regions. The 
Pesultant velocity potential distribution is used to deter- 
mine the pressure loading over the wing through the rela- 


maonship 


eee= =cp.,\U = + iw)? ee) 


co 


taken for the harmonic case from equation (2.18). 
meeviake Eifect 

One of the historic drawbacks to the potential 
approach to lifting surface theory has been the necessity 
of including the wake region in the solution of the integral 
equation (3.1). This equally complicates the discrete 
element approach since the matrix equation -(4.2) is required 
memetcluce the effect of the wake potential dipole grid. 
Even though the effective wake is terminated at some repre- 
Pewraci ve length behind the wing, as is commonly done in 
steady flow problems, the size of the resulting potential 
strength vector and kernel function matrix would severely 
Beoeeercb, Or entirely preclude, the use of this method even 
on large modern computers. 

Porlewing a2 Sugeestion by Professor R. E. Ball of the 
Naval Postgraduate School, a careful examination of the 
boundary conditions governing harmonic wing motion, as 
depicted in Figure 3, was made. In the wake region the 


pressure is zero, but the velocity potential is not zero 


oe 
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Mecause of the discontinuity in the u component of the per- 


meioavion velocity. $ and p are related by 
A Q 


In the wake, therefore, this relationship becomes 


U, <2 + iwd = 0 (4.5) 
Ong 

od _ WW 

TT ig 


which can be integrated to give 
Ww 
=, 


UL (Xx -X) 

Do ae. © C426) 
a ime menserengeth of the wake velocity potential at 
(xX 3¥4 20); and ¢ is the strength of the VieMlocieymmOtemt iad 
at the boundary between the wing and the wake (X3¥5 20) 
a@€long the line of integration Vavae In the discrete poten- 
tial element formulation, 2m becomes the strength of the 
last wing dipole on the appropriate chord line, with Ce) 
the distance between this dipole and the wake dipole 9, 
as is shown schematically in Figure 4. 

The wake potential strength distribution is therefore 

Bepressible as a function of the wing potential distribution 


fwiemcan be included in the wing doublet kernel functions. 


This relationship can be represented by 


ere Ullal¢ } 
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where {o} is the wing velocity potential vector and the 
appropriate elements of [A} express the relationship of 
equation (4.6) with respect to the wake potential vector 


{to}. The downwash matrix equation (4.2) therefore becomes 


{w} 


[K,] {6} + [K,] {4} 


[K,] {o} + [K,] [A] {o} 


[K, +K,A] to} 


mifewstZze Of this matrix equation has thus been limited to 
meee necessary to represent only the wing grid, making the 
Giscrete potential element approach much more tractable on 
Che computer. 

pee liatrix Equation for Symmetric Motion 

In considering a wing undergoing harmonic oscilla- 

tions which are symmetric with respect to the x-z plane, the 
Wing/wake planform can be represented as shown in Figure 5. 
The x-z plane, being a plane of symmetry, represents a no- 
flow surface, or reflection plane, for the disturbances 
feeseea OY the symmetric motion of the full-span wing. View- 
ing the physical problem in a different, but equally valid 
light, the half-span wing @, with its root along the x axis, 
can be considered undergoing oscillations in the presence of 
an infinite wall coincident with the x-z plane. Therefore, 
the motion of the wing need only be prescribed for wing Q, 
while the function of the virtual wing @ and wake @is to 
establish the no-flow condition at the x-z plane. 

The full wing/wake system is represented by the matrix 


equation 
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a 


Wy Ky, | Kyo | Ky3 | Ky >> 


a0 | See Sacer eae 
| >y 


where subscripts denote the areas of Figure 5 and Ks is the 


kernel function matrix for downwash at points in areaQ@ due 


femoeublets located in area(j). For symmetric motion 
(@,} = {6,} and {64} = (65) 


so the required downwash equation is reduced to 


oT 


- 1 
{tw} = [Kj ,+K, 3 ! Ky otk, 4] 


In the previous section, the relationship between wing and 
wake potential distributions was developed, and can here be 


expressed as 
{o,} = [Al {6} 


Therefore, the final wing downwash matrix equation has the 


mem Of Equation (4.2) 


{w, } = [K] {$5} CH) 


but the kernel function matrix is here formed from the 


following expression 


Thus, through use of motion symmetry and the boundary condi- 


tion in the wake, the effective integration of equation (4.7) 
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meed only be carried out over the half-span wing surface Q. 
Acknowledgement of the full physical problem is obtained 
eeeoueh formation of the kernel function matrix via equation 
(4.8), which incorporates effects of the potential dipole 
memos 19 all four areas of the wing and wake depicted in 
Figure 5. 

It should be noted that antisymmetric motion of the wing 


about the x-z plane requires that 


(¢,} = -{6,} and {9} = -{6,) 


The final kernel function matrix is therefore formed by 


[K] = (K,,] - (K,] + (K,5-K), 104] 


Since any general harmonic motion of the wing can be 
expressed as a combination of symmetric and antisymmetric 
modes, this approach has general application and need not 
be restricted to the symmetric case considered in detail 
here. 
4. Doublet Singularity 

The location of the collocation points required to 
wmlety the no-tlow condition in lifting surface theory is 
historically based on two-dimensional steady flow theory, 
Vortex lattice methods such as Rodden's [38] arrange the 
Membre x Strip on the local quarter chord of the control box, 
with the downwash collocation point centered on the local 
three-quarter chord line as in two-dimensional thin airfoil 
maeory. Houbolt [6] ee ese Woti@ebiits Local quarter pnoed., 


Bmiree—-Guarter ehord control box grid in conjunction with 
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concentrated pressure loads. This type of grid network was 
tried by the author in an early form of the wing analysis 
memoucer program with unsatisfactory results. The approach 
Mmaaetound not to converge, but to be very sensitive to grid 
mmc chat iS, to the distance between the potential dipole 
memos associated control point. This is much like the 
sensitivity that a continuous loading method experiences 
Meem a control point is located too close to a wing edge 
orl). 

Meee onleral or collocation points were subsequently 
placed at the center of the wing control boxes (Figure 2) 
@emmeldent with the potential doublet locations. In this 
Way, the above mentioned sensitivity to grid spacing was 
removed, but the value of the upwash of a doublet at its 
Own control point had to be determined. As can be seen 
from equation (4.3), the kernel function (Ky) is singular 
meee Goublet location (i=j). This singularity in the 
upwash from a doublet is pictured in Figure 6a, where the 
feuplet produces an infinite upwash at its location, but 
finite and decreasing downwash in the plane of the wing as 
the distance from the doublet is increased. 

If the dipole is considered within the framework of the 
discrete element grid where dipole strengths and downwash 
velocities are held constant, or averaged, over each control 
box, a cross-section of the dipole flow pattern would appear 
as in Figure 6b, where the infinite upwash at the doublet 


has been replaced by a finite value Wo: To determine a 
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value for Wo Woikenm adequately represents the singularity 
Berength within the framework of the discrete element 
meeroximation, the law of continuity was employed. The 
meer tid Outflow in the discrete representation is WoAo> 
where Ay is the area of the doublet control box. Modeling 
the entire x-y plane with control boxes without regard to 
Manie/wake geometry, the total amount of fluid passing back 
through the x-y plane is Z w A,» where w, is the average 


feroctuy over the control box with area A: Therefore, 


Bemoinuity requires 


and if the A's are chosen so that A, = An: the upwash 
velocity is determined by 
Wo - a we (4.9) 
In actual practice, the summation of the downwash velocities 
through the x-y plane is necessary only to some finite radius 
from the doublet, since the inverse proportionality of the 
Peanet function with distance from the doublet causes the 
value of We tO converge within a reasonable summation. 
peepeccuron and Wing Coefficients 
To obtain section lift and moment coefficients, the 
chordwise pressure distribution at each spanwise station 
was integrated in the manner employed in thin airfoil theory. 
First a coordinate transformation of the chordwise variable 


is made as follows 
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: | e~ \ 
=o. fe 
O 0.5 1.0 
(L.E:) Cee 


Xf. = 5(\-cos 6) 


Meemmorcosure coefficient is then expressed in terms of a 


Fourier expansion of the new variable 


"om 8 


fs) 
c (0) = ec ar 
* ) 5 cot 5 


ec sin né (4.9) 
O n p 


il n 


where the first term accounts for the leading edge singu- 
larity, while the trailing edge slope singularity is 
acknowledged by the infinite series. The section lift is 


obtained from 


C, a2 He 0) 





Substituting the Fourier expansion for C,, B5qG| joo i iee@unabayes Pi elavs: 
coordinate transformation, equation (4.10) becomes 
e E 
ea ee 00 Ph oT 
Cy = awe Girces oda + 2 = [ecimaene siniedg 
oi (4,12) 


OtewcG the orthogonality of the sine function, this integra- 


tion yields 


= 5 (ee + Sc ) Cae 


The section moment about the leading edge is defined by 
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Cc =-f ee 
ae QO 


d(=) (4.13) 


OQ |x 


Dp 


Performing the same transformation and integration as for 


the section lift coefficient, the following results 


is 1 
em ~~ B Spt Sp” 2 Sp.) (4.14) 


meamsterring this moment coefficient to the mid-chord point 


requires the further calculation 


The integration of section lift and moment valuesto wing 
@eeiticients is performed in a similar manner. Here the 
coordinate transformation is made to the spanwise variable, 


€é 


euem that 


y/b = - (l-cos 6) 


Mae product of the section chord and the section lift 


CeGefficient is then expressed as the Fourier series. 


Q oi 
ce, = c, cos ~ + ££ ec, sin née Ce, alicia 
NE Mee 2 Areal L 
Pemrorming the integration 
Ce : ec, d(¥) (4.17) 


of the Fourier expansion of the section lift in the trans- 


formed coordinate system, the following results 


_ 8 2 TT 
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The same relationship is obtained for the wing moment 
coefficient when the sectional moment is expressed as in 
Mee Pourier expansion of equation (4.16). For swept wings, 
meeesectional moments were transfered to an axis extending 
meem tne half—-root-chord point, prior to developing the 
Pourier series. 

The infinite series of equations (4.9) and (4.16) are 
of course terminated at the number of chordwise and span- 
wise points respectively of the wing discrete element grid. 
This places a minimum limit on the number of chordwise 
and spanwise control points which are required to achieve 
valid sectional and wing coefficient values, a factor which 


will be discussed in Section VIL. 


B. WING/BODY ANALYSIS : 

ieee sasite Formulation 

Inclusion of a finite radius body in the nonsteady 

lifting surface analysis is equivalent to adding one more 
boundary condition to the problem formulation: the require- 
fier tor no flow Chrough the body surface. For the analysis 
Pemsued here, the body will be stationary with regard to 
the perturbation motion, and will be idealized as an 
infinitely long cylindrical surface with axis coincident 
with the undisturbed flow. The cantilevered, midmounted 
Woe is harmonically oscillating within the limitations of 
small disturbance theory previously developed. 

In steady flow amare tee tier Dectwwenrhocts have been 


traditionally handled by singularities, matching the wing 
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meaeularity distribution, at image points within the body in 
the wing plane [39]. This method will satisfy the nonsteady 
boundary conditions only in a quasi-steady sense, because 
Mimpinase differences between disturbances at the body sur- 
face caused by a wing singularity and its corresponding 
image. Therefore, to model the body in the nonsteady 
meeorem, a2 system of singularities are placed on the body 
feaetaee establishing a grid similar to that for the wing. 
These curvilinear panels each have an harmonically oscillat- 
ing singularity at its center, coincident with a normalwash 
control point. Thus, in effect the wing grid is merely 
extended over the body surface as shown in Figure 7. 
imieorder LO handle the more complex geometry of the 
wing/body problem, a coordinate transformation from the 
rectangular (x,y,z) system of the wing analysis is made to 
a cylindrical (r,68,x) system, where the undisturbed flow 


Girection (x axis) is common. Thus the transformation is 


defined by 
x = xX 
y =r cos @ (4.19) 
7e@—- © Sin 6 


Figure 8 shows the wing/body configuration with the boundary 
Semeitions in the cylindrical coordinate system. The down- 
wash velocity over the wing is now related to the velocity 
Pemcttvial distribution by 


cx) (4.20) 
a S=ONoOr TT 


= 
cD 

ll 
a= 
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WING/BODY CONFIGURATION IN 
CYLINDRICAL COORDINATES 





and the no-flow condition for the body surface is stated as 


u,, == oa =o) Cle 2) 
r=r 
B 
where Pp is the body radius. As in the wing analysis, the 


wing downwash velocity Ug is determined from the type of 
surface motion prescribed for the problem. 
eee Governing Matrix Fquation 
The governing matrix equation for the wing/body 
problem takes on a much more formidable appearance than the 


basic equation (4.2) for the wing along 


Ye) (Sy | Swel |?w 
—-t= |—--—-4--] {-- Gie2)) 
Ur! [Kew y SB] [%, 
where {ug} = wing downwash vector 
tu} = body normalwash vector 
{oy 5 = wing singularity potential vector 
{dp} = body singularity potential vector 


The kernel function matrices are defined as follows 


Ky = downwash on wing due to wing singularities 
Kup = downwash on wing due to body singularities 
Kew = normalwash on body due to wing singularities 
Ke = normalwash on body due to body singularities 


rece Solution of equation (4.22) for the velocity potential 
distribution would be, limited due to the computer storage 


requirements of the large complex kernel function matrix. 
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However this approach is not feasible at any rate, since 
meeererne! function matrix is ill conditioned in this form 
mmmerocs NOt lend itself to efficient inversion. 

To attack the problem, equation (4.22) is separated into 


two matrix equations 


(oy Ee} + IK plop} 


feu} 


. [Kylie ,t + [Kp] 1p} 


The second of these can he solved for the body potential 


sori bution 


{og} = (KpI7> [u,-K, $y] (4.23) 


meee 1S then substituted into the first equation. 
MPS Ke Me) PGK (Kalo tu.} = (K,) (Kal 1 {Kayl (oy) 
6 W W WB B r WB B BW W 


This may be rewritten in the following form 


ai _ Sal 
K ue 2 TK a Sele Kpy! toys Cie) 


{ug-Kypkp 


The left hand matrix is a modified wing downwash vector 

incorporating body effects, as does the modified kernel 

function matrix on the right hand side of equation (4.24). 
In the analysis considered here the body is steady at 


Mero angle of attack to the main stream, while the wing is 


midergsoing nonsteady motion, so that 


{u,} = {0} 


Therefore, equation (4.24) becomes 


Bi 





ills II Ky-®woke Koy] 25) 
which is in the same form as equation (4.2), but in which 
the modified kernel function matrix incorporates the body 
Beuncery condition. 

3. Symmetry 

The individual kernel function matrices of equation 
(4.25) incorporating body effects can be quite large when 
Semvrol points are placed around the entire circumference 
Geeeene body, severely limiting the utility of this method. 
Pemever toese Matrices may be reduced appreciably in size 
through use of symmetry. Considering the representation of 
Figure 9a, it can be seen that for wing motion symmetric 


about the x-zZ plane 


4 


Cty} = (Oy 35 Cop} = Cop ts Lop} = op 3 (4.26) 


where {$,, | stands for the doublet strength distribution on 
the ec ctive wing, and {o, I SOendemt orm thessineulari by 

eemeeneth distribution on Hae body in the respective quadrant. 
Jt can also be seen that due to the antisymmetric nature of 


the wing doublet flow with respect to the x-y plane 


(¢p) = - pts (op d= - (4g (4.27) 


Since the purpose of the body singularities is to counter 
the flow of the wing doublets through the body surface. 
Considering first the body kernel function matrix (Ky); 


which is defined by the relation 


Sis 





{u_} = [K,]{¢,} 


or 


Lu 


Jim 


3 3 


symmetry allows this formulation to be reduced to 


} Cis) 


ues Cn eee = K : 


}i¢ 
B, By, “Bz By °'B 


In the case of the wing-body kernel function matrix (Kip), 


this same consideration produces 


{ug} = [ik K -K -K 119, } GimeD 


+ 
WB, WB, vee WB), 1 


Control points need only be placed, therefore, over the 
first quadrant of the body surface, reducing the size of 
the body effect kernel function matrices of equation (4.25) 
by one-fourth. The surface over which the flow conditions 
are specifically satisfied are indicated by the solid line 
of Figure 9b, which then satisfies the boundary conditions 
on the whole wing/body surface through symmetry. The body- 
wing kernel function matrix (Kaw) follows the same develop- 
ment as in Section IV.A.3 for the wing alone and has the 
Same form as K,, given by equation CH tay 

It should also be noted that for antisymmetric wing 


motion the kernel functions of equations (4.28) and (4.29) 


would have the form 


l= AK, -- K 


a + K 2 ~ Ky] 


BM, 


Pee + (Ke Jig, f+ [Ko J{¢, } + [Ke ]{¢ 
B, ''B, B. B. B B By By 
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Since any general harmonic motion of the wing can be 
expressed as a combination of symmetric and antisymmetric 
motion, this method, as in the wing alone case, need not 
be restricted to the symmetric motion case considered in 
detail here. 
4, Effective Body Length 

The body considered in this analysis is idealized 
as an infinitely long cylinder with axis aligned with the 
undisturbed free stream flow. This is not an unrealistic 
imi'mcrcetion, Since in most practical cases, the center part 
of a fuselage is nearly cylindrical and the fineness ratio 
of the fuselage is large enough so that the flow at the 
center part is nearly the same as for an infinitely long 
cylindrical body. Experiments in steady flow have proved 
that beyond an effective body length, the lift distribution 
of the wing/body combination is independent of the body 
Memeci [40]. Steady flow analyses, such as Woodward's [41], 


employ a "wing-body interference region," where the body 
no-flow boundary condition is explicitly applied a finite 
distance upstream and downstream from the wing root section. 
A similar effective body length is modeled in this 

approach, where the body singularity grid extends from a 
finite distance upstream of the wing root to a point down- 
stream of the effective wake. The actual body length which 
must be modeled in order to include all interference effects 


18 discussed in Section VI. Since the kernel function 


velocities decrease approximately as the distance from the 
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mmmeularity cubed, it is not unreasonable to assume that 
interference effects are concentrated close to the wing/wake 
area. 

pee DbOOy singularities 

The acoustic singularity used to model the wing was 
a potential doublet, due to the requirement for a pressure 
Gitferential across the wing surface. No such requirement 
feerous tor the body singularities, so that an harmonically 
oscillating source could be employed as well as the dipole. 
However, in the development of the wing computer progran, 
meocecame apparent that the solution to the matrix equation 
(4.4) was very sensitive to the value placed on the kernel 
function upwash singularity. Approximate methods used to 
obtain an average value of the upwash did not produce 
acceptable loading distributions over the wing. Only when 
this upwash was determined numerically, as indicated in 
Section IV.A.4 were good results obtained. Unfortunately 
a source singularity, while having a somewhat simpler form 
of kernel function, can not be analyzed by this same type 
of numerical development. 

The body was therefore modeled by a network of harmon- 
ically oscillating potential doublets, as on the wing, with 
meceOraenvucd in the radial direction. Each doublet is at 
the center of its control panel, coincident with a normal- 
Wash control point. Determination of the normalwash (radial 
velocity 2 at the doublet location is obtained, as with 


the wing singularity, from a consideration of continuity. 
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Mie coOval fluid outflow from the doublet is U, Ay where 
0 
Ay is the area of the doublet control panel on the body. 
The amount of fluid passing back through the body at a 
control panel away from the doublet is (u_) A where n 
Ponm nm 
indicates the axial location and m the circumferential loca- 


memermeor this control panel. Continuity is then expressed 


fOr this discrete element representation as 


00 M 


u,, a + fF % (u,,) oe = 0 
O n=1 m=l1 nm 


where M equals the number of panels located around the 


circumference of the body. If the Bee ee meMesenins Omelans 


am = Ay the value of the doublet singularity is given by 
oo M 
UL ee 2 (u_.) ‘ (4.30) 
O n=1 m=1 nim 


ieee cOnvINnuity equation states, in effect, that all the 
Mmueomwoilehny leaves Che doublet in the radial direction must 
Pess back through the body surface prior to returning to 
the doublet. As in the wing singularity case, the summation 
Seenormalwash in the axial direction can be terminated after 
a reasonable distance in both the upstream and downstream 
sere c tions. 
6. Kernel Functions 

The elements of each of the kernel function matrices 
Pemined in Section iV¥.8.2 represent the normalwash on either 
Wing or body surface due to a harmonically oscillating 
doublet of unit strength which satisfies the governing 


linearized differential equation (2.12). Doublets on the 
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wing are of course oriented in the 8 (or z) direction, 
while those on the body have their axes in the radial 
meee ction. 

The formulation of these kernel functions is developed 
in Appendix A. They are summarized below with control 


a ee ine 


Meno at (r,8@,x) and doublet located at (ro 38); = 


meay radius is r 








B° 
Ga) K, - Downwash on wing due to unit wing/wake singu- 
larity. 
- Be fe 100) W 
Ug = az (1 +i —S5 Roexp {iB [Mi (x-x,)-R]} (4.31) 
UaR a8 a8 
where 


Po 2 ap 
R ay eran) te (te Wiens Aire) ; 








ei) K,,, - Normalwash on body due to unit wing/wake 
Singularity. 
2 Bye 
= = sclial [8 Qo SS) (14+4 —5 R) 
4aR R a_B 


Loo) 





2 
W JR W 
+ (>) [rpR ge] exp {1 ee: [m.(x-x,)-R]} (4.32) 


co 


where 
Z 2 a a a 
mS NOR) + 8 (rptr. a 2r orp cos 6) 


B + Po cos 86) 
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In both the above formulations, the minus sign represents a 
singularity on the 6=0 wing/wake, and the plus sign a 
singularity on the 6=n virtual wing/wake. 
(Gi ) Kp - Normalwash on body due to unit body singu- 
penta Wa 


O 
_ =f ie! 7 OR 
u,, = aoe | eosce-8,) + = rp (1 - cos<6 87) a 








. Ww eee oR 
+4 5 - ay Rr, (1 - cos<8-6 >) ae | 
a8 a_B 


co 





exp [4 3 [m,,(x-x,) - R]} (4.33) 


where 


R= J (x-x5)° + 28° re (ES COSSISE 2) 


o 
oho 6 
ae — rp (1 - cos<6-6 >) 


(iv) K,, - Downwash on wing due to unit body singu- 
lene es 


_g° 


ie ’ 3 oR 
a r sin 6 + = (r,- r cos 8_) | 








2 
a+ 3 R - (—) R( rp - r cos 8) | 
a_B elas 


co 





exp {4 a [I (x-x,) - R] (4.34) 


_ Oo Ae Be 
R= V (x-x,) mp (re een cos 6) 


aR Be ea 
20 7 9: 5): aie 


where 
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V. DESCRIPTION OF COMPUTER PROGRAMS 


ms GENERAL DESCRIPTION 

Listings of the two computer programs developed and used 
in this investigation are reproduced in Appendices D and E. 
The programs are written in standard FORTRAN IV language. 
Calculations were performed on the I.B.M. 360/67 computer 
at the W. R. Church Computer Center of the Naval Postgraduate 
Bemool. Object codes were obtained with the I.B.M. G-level 
eomoller. 

Both the wing and the wing/body programs are arranged in 
the same general format. The MAIN program reads the input 
data and establishes the wing or wing/body control Oueele 
Senor more subroutines are called which establish the 
kernel function matrix elements. MAIN calculates the wing 
Gdownwash velocity vector and calls subroutine COMAT to 
solve the matrix downwash equation (4.2) for the wing velo- 
meee pOLeHoial Gistribution. MAIN finally calls subroutine 
PRES which calculates the perturbation pressure distribution 
over the wing by applying finite difference approximations 
GO equation (4.5). PRES then calls subroutines SECLM and 
WINGLM to integrate this distribution to obtain sectional 
and wing forces and moments. These results are printed by 
PRES and control is returned to MAIN for the reading of 


data cards for the next problem. 
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B. WING PROGRAM 

The wing program computes potential and pressure 
loadings due to harmonic motion on general planar wing 
configurations from rectangular to arbitrary sweep angles 
Tor both the leading and trailing edges. The program is 
maemeeleved, however, to constant Sweep angles, and to 
a finite tip chord with minimum taper ratio of about one- 
meee o. Lhus delta wing configurations are excluded. This 
restriction is caused by the method used to integrate The 
wmerowlse pressure distribution to obtain sectional lift 
and pitching moment, discussed in Section IV.A.5. The 
mmeeram provides for a maximum of 100 control points on 
Pieswing and ten control points in either the chordwise 
or spanwise directions. Storage requirements are the 
equivalent of 18,000 single precision complex words, or 
144,000 bytes on the I.B.M. 360, with a maximum run time 
of approximately twelve minutes for one problem. Provision 
mmade 1Or the running of successive problems for as many 
sets of input data cards as are provided. 

Figure 10 is a diagram of the wing program, while sub- 
routine flow diagrams are presented in Appendix C. Sub- 
routine DINCO forms the kernel function matrix from 
equation (4.8), with the matrix elements defined by equation 
22) The COMAT subroutine solves the linear complex 
matrix equation (4.2) by the Gauss-Jordan method with total 
pamiouing. This subroutine was written for systems of real 
Pemioetons by Mrs. Sharon Good, David Taylor Model Basin, 


as published in Ref. 42, and modified to include the complex 
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@epabtlity by Mr. Hellmut Golde, Department of Electrical 
Engineering, University of Washington. COMAT was further 
Mear171ed by the author for particular application to the 
wing and wing/body programs. 

Mmieur instructions are presented in Appendix B. Any 
number of spanwise and chordwise control points, up to the 
maximum of ten each, may be specified for a semi-span wing. 
Mhe wing is modeled by an equal number of chordwise control 
points at each spanwise station. The program is limited 
by theory to the subsonic flow regime, within which any 
mode, amplitude, or frequency of harmonic wing motion may 
be specified, including the steady case. Results printed 
by the program for each spanwise station are: 

G) Potential doublet strength distribution 

Gia) Pressure coefficient distribution 

(iii) Section lift and pitching moment coefficients 
Maeecditi0n, wing lift and pitching moment coefficients are 


presented. 


C. WING/BODY PROGRAM 

The wing/body program incorporates the body surface 
boundary condition into the harmonic motion analysis of a 
rectangular wing planform. As in the wing program, a maxi- 
mum of 100 control points may be used to model the wing, 
while up to 130 points are allowed inthe body control panel 
network. Storage requirements are the equivalent of 43,500 


single precision complex words, or 348,000 bytes on the 
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I.B.M. 360. A maximum run time of about 20 minutes for 
the wing/body problem is required, while a wing only solu- 
tion (body radius equals zero) is obtained in less than 
seven minutes. As in the wing program, successive problems 
fermediitiferent wing, body, and flow geometries may be run. 
Figure 11 diagrams the wing/body program, with indivi- 
dual subroutine flow charts presented in Appendix C. The 
@ammeenece coecificient matrices defined in Section IV.B.2 


meen tormed in Che following subroutines: 


i) DINCO - Dy from equation Ge) 
Git) Vines = Kwp from equation (4.34) 
(i137) URAD1 - Kaw from equation (4.32) 


inv ) URAD2 - Ke from equation (4.33) 


The matrix kK. is inverted by subroutine COMAT and the modi- 
fied kernel function matrix of equation (4.25) is formed 
in MAIN. COMAT is again called to solve the resulting matrix 
Eduavion and PRES performs the same function with the same 
Mew printouts as in the wing program. Subroutine DINCO, 
while performing the same function as in the wing program, 
is here in a more simplified form because of the restricted 
wing geometry of the wing/body program. This accounts for 
@eae reduced running time of wing only problems in this 
program as compared to the general wing program. 

The wing/body program uses the same wing grid and flow 
geometry parameters as the wing program. Any number of 
meay control points may be specified up to the maximum of 


130. An effective body length of one wing chord upstream 
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FIGURE 11 
WING/BODY PROGRAM DIAGRAM 
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from the wing leading edge to one wing chord downstream 
from the termination of the effective wake is modeled by 
meee program. input data instructions are presented in 


Appendix B. 
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VI. RESULTS AND DISCUSSION 


A. WING ANALYSIS 
meecomparison with Lifting Surface Theory 

In order to investigate the validity of the discrete 
DOvential element approach to lifting surface theory, the 
wing program was run for a wide variety of wing/flow geometries 
and types of oscillatory motion. Examples were picked which 
could be checked against previous work reported in the 
literature, representing a variety of lifting surface theory 
approaches, as well as the relatively limited amount of 
experimental results which have been obtained in the unsteady 
field. 

The figures presented in this section are a representative 

mammary Of this investigation. The subsonic flow regime was 
Bovered irom the incompressible M_=0 to the high subsonic 
ieee ?, While the range of frequencies wricd fron the steady 
case to the relatively high reduced frequency of k=1.2. 
The wing planforms considered had a range of aspect ratios 
from one to six, sweep angles from zero to 45 degrees, and 
taper ratios from one to one-half. The types of symmetric 
harmonic motion considered were: 

Coe trching = rigid body oscillations of the wing 
about a spanwise axis perpendicular to the flow direction; 

(11) Bending - oscillations of the wing as a beam canti- 


levered at the root chord in the first natural mode of bending. 





mao) Plunging — rigid body oscillations of the wing in 
the z direction at a mean angle of attack of zero; 

cy) Mentored body roll oscillation of the wing 
mepeey Supported from a chordwise axis inboard of the root. 
steady data were obtained with the wing analyzed as a rigid 
body at a fixed angle of attack. The results have been pre- 
sented as chordwise pressure distributions, spanwise lift 
Pm@omoleoCching moment distributions, and wing coefficients 
Memeroed With respect to reduced frequency. 

In each of the figures, the coefficient amplitudes are 
normalized with respect to the motion as follows: 

(i) Pitching - angle of attack amplitude; 

ci) Bending - wing tip angle of attack amplitude; 

(iii) Plunging - vertical motion amplitude; 

(iv) Flapping - amplitude of Pesos angle. 
these represent the conventions used in the applicable 
references, however the phase angle relationship were con- 
verted to the coordinate convention used in eae paper where 
Pamererences occurred. 

Ser releavion Of the wing program results with existing 
Peesine Surface theory was in general quite good. Devia- 
tions appeared to come from the different methods employed 
momnmandle the wing edge singularities in the pressure 
Gistribution. The lifting surface approaches assume the 
horm Of these singularities a priori in the choice of their 
Bencane functions, as discussed in Section III.B. The 


discrete element approach, on the other hand, assumes no 
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loading profile, but obtains proper potential and pressure 
Meeeributions through inclusion of the boundary conditions 
mee ne problem formulation. As discussed in Ref. 1, much 
eitfort has been devoted in lifting surface theory to 
developing the most numerically efficient and accurate 
Poading functions, as well as to improve the numerical 
methods of handling the singular kernel functions. 

The first four figures compare wing program results with 
one of the most recent lifting surface theories. This 
advanced kernel function method, based on the acceleration 
potential, was presented by Laschka in 1963 [5] and further 
developed by Laschka and Schmid for interfering planar 
surfaces in a 1967 paper [43]. Figure 12 compares chordwise 
Pressure distributions on a swept, tapered wing, undergoing 
pitching oscillations in low subsonic flow with Laschka's 
results [46]. Good correlation is evident at each spanwise 
station, with the only variance being a small difference in 
the shape of the pressure distribution near the leading 
edge, as previously discussed. 

Figures 13 and 14 present spanwise loading on a 45 degree 
swept constant chord wing in both the pitching and plunging 
modes. Results are compared with Lashka's work [5] for both 
ten subsonic (M = 0) and M. = 0.8 Plows correlation 1s 
again good. The discrete potential element approach appears 
to overestimate the amplitude of the lift and pitching 
moment slightly, especially near the wing tip. In fact this 


latter trend appears in almost all the results of the present 
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method. The explanation again would appear to lie in the 
meowood Of handling the wing tip pressure slope SineuLaracy. 
Weeety 15 included implicitly in the boundary conditions. 
Further development of the discrete potential element 
meproach would definitely have to include investigation of 
the adequate recognition of this wing tip condition in the 
mmoolem formulation. 

Hn Ref. 3/ LaschKa compares his results, for the wing 
configurations of Figures 13 and 14, with the work of Pao 
Tan Hsu [44]. It should be noted that the wing program 
results agree more closely with Laschka's data, while fall- 
ing between Laschka's and Hsu's results except for amplitudes 
at the wing tip. The theory here is further compared with 
experimental results obtained by Laidlaw [1B with fair 
morrelation. 

Figure 15 compares wing program spanwise loading with 
Laschka's results from Ref. 37 for a steady swept, tapered 
[meetin low subsonic flow. The lift and pitching moment 
coefficients are presented on an expanded scale, with the 
Same type of comparisons already noted. 

The results presented in Figures 16 through 18 are note- 
worthy for the Mach number range considered (M,=0.24 to 0.9) 
and for the consistent experimental data presented from the 
work of Lessing, Troutman, and Menees [47]. Spanwise and 
chordwise loadings are plotted for an aspect ratio three 
rectangular wing in the bending mode. Theoretical .comparison 


is made to the lifting surface results developed by Lessing 
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meal in 1960 from the original kernel function methods of 
Refs. 25 and 48. 

At a Mach number of 0.24, the chordwise pressure loadings 
Meearmed from the two theories are almost identical. The 
spanwise lift and moment curves show a somewhat heavier 
loading concentration near the wing tip for the discrete 
potential element method, as in the preceding figures. This 
trend also appears at Mach numbers of 0.7 and 0.9. The 
Beano. ( Chordwise pressure distributions show a difference 
in the representations of the leading edge pressure distri- 
butions outboard of n = 0.5. These deviations can probably 
again be laid to the different methods employed to handle 
Che wing pressure singularities. However, the wing program 
results agree more closely to the experimen vat data of Ref. 
47, than does the kernel function approach. No explanation 
can be found for the deviation of the pressure phase angles 
in the wing trailing edge area for both the 0.24 and 0.7 
Mach number cases. It should be noted, however, that the 
experimental points again correlate more closely with the 
meme program results in this area. In fact, the close cor- 
relation of these experimental data with the results of the 
present method is quite gratifying. 

It is interesting to observe that a shock wave may have 
started to form on the wing in the 0.7 Mach number flow down- 
stream of the 60% chord line, as indicated by the drastic jump 
in the pressure phase angle measurements at this location in 
Figure 17b. The theoretical analyses would not, of course, 


reflect the existence or the shock wave. 
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The closeness of the two theoretical approaches at 
foo? 25 also noteworthy. As discussed in Section II.C, 
Meenlinearization of the basic problem for perturbation 
analysis would appear somewhat questionable this close to 
the speed of sound. However, the consistency of the results 
would seem to indicate that the linearized subsonic theory 
Meesctill valid at this high subsonic Mach number. 

mieeure WO analyzes the pitching motion of an aspect 
ratio two rectangular wing. Chordwise and spanwise loadings 
are compared with the experimental and theoretical results 
of Laidlaw [45]. The latter represents a numerical treat- 
ment of the rectangular wing aspect ratio theory of Reissner 
[49, 50] developed in 1947. Considering the fact that 
Reissner's approach basically involves applying Correction 
imaerors tO two-dimensional results in order to account for 
finite aspect ratio, this theory agrees quite well with The 
discrete potential element approach. The experimental data, 
although not as consistant as that of Lessing [47], agrees 
reasonably well with the Wing Program results. Figure 20 
presents spanwise loading for the same wing in plunging 
ieemron. Correlation of pressure distributions between 
theories and experimental data for this mode of nonsteady 
MmoG1on is the same as for the pitching case. 

The frequency response of a rectangular wing in the flap- 
ping mode is presented in Figure 21. Comparison is made to 
experimental and theoretical data by Woolston, Clevenson, 


and Leadbetter [51], who employed a basic kernel function — 
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meemoach. Correlation of the wing program results with the 
experimental data is reasonably good, although there is some 
difference in the phase angle values. However, the varia- 
tions of both lift and pitching moment with reduced fre- 
quency compares very well. Woolston's theoretical points 
coincide essentially with the wing program curves except 

for the pitching moment phase angle. No explanation can be 
Mound for this difference. Laschka [5], in comparing his 
results with those of Woolston, produced virtually the same 


curves as the discrete potential element program. 
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Eee 6 COnvergence 

Figures 22 through @ present spanwise loading for 
meeee Gitferent wing planforms as a function of the singu- 
darity grid density in order to show convergence of the wing 
program. Three types of motion are considered, with wing 
memos numbering from 25 to 100 points. Figure 25 shows 
convergence of the wing program for the steady case on an 
aspect ratio six rectangular wing. This wing planform is 
used in the basic reference for the wing/body analysis in 
steady flow. The ordinate of these graphs has been greatly 
expanded to permit a good evaluation of convergence. 

All of the examples are seen to converge at about the 
same rate. At 50 control points on the wing the maximum 
deviation is less than four per cent. At 60 or more points 
the results are virtually coincident. The shape of the 
Pemerol bOx would appear not to be a factor, since varying 
the grids on the individual wings of Figures 22 through 25 
also varies the shape of the individual boxes. However, it 
would seem prudent to maintain the control boxes at reason- 
able aspect ratios (less than about ten) to assure valid 
results. Another restriction on the wing grid size is the 
requirement to integrate the pressure loading to obtain 
section and wing coefficients. A minimum of six points in 
either the chordwise or spanwise directions was found neces- 
sary to adequately define the loading distribution for the 


integration subroutines. However, the basic requirement for 
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60 or more points in the wing grid makes this latter 
meee lcti0on a factor only for higher aspect ratio wings. 

A second type of restriction on the wing grid size is 
caused by the accuracy of the inversion method of subroutine 
COMAT. Solution of the downwash matrix equation is accom- 
plished by the Gauss-Jordan method employing total pivoting. 
With the relative magnitudes of the kernel function matrix 
Pwemenbs Occurring in this type of analysis, the inversion 
fee@ecess Starts to lose accuracy with a system of equations 
or order greater than about 110. Computational accuracy 
feeeeveori tied by Substituting the velocity potential vector 
solution into the downwash matrix equation, and obtaining 
Emoownwash vector to compare with the original input. Checks 
Peeeees SOlUuULiOn for grid densities up to 150 points were 
Meade, 10r various wing/flow conditions, with consistent loss 
BPeeoccuracy. Thus, the program has been limited to wing 
meteor 100 control points or less to ensure accuracy. 

This also limits the aspect ratio of the wings that may be 
analyzed; however, for the purposes of this report, aspect 
ratios of six or less were well within the capability of 
the program. If added capability were required, a more 


sophisticated inversion routine would have to be developed. 
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3. Wake Effect 

Figures 26 and 27 summarize the effect of the finite 
wake length considered in the singularity grid network as 
@mscussed in Section IV.A.2. Seven configurations of wing/ 
flow geometries were analyzed, covering the range of condi- 
tions considered in this paper. Wake effects were included 
Pena Maximum distance of ten root chord lengths downstream 
from the wing trailing edge. In all cases, wing loading 
had converged to within one per cent of amplitude and one- 
half degree of phase angle by four root chord lengths of 
effective wake. This correlates with Haviland's work [8], 
in which he reported that the effect of the wake on the wing, 
for a rectangular planform in the steady case, was determined 
meomwercohnin one per cent in five chord lengths. 

Figure 26 presents a specific case of wing coefficient 
variation with effective wake length. The maximum deviation 
of the wing coefficients, with respect to an effective wake 
meee Ot four root chords, is shown in Figure 2/7. From 
the latter figure, it can be seen that the near wake (within 
one root chord length of the trailing edge) has a very great 
effect on the wing loading. The wake influence then dies 
out quite rapidly as the effective termination point is 
moved downstream. Therefore, the finite wake, taken into 
account in determining the aerodynamic loading on the wing, 
need only extend a realtively short distance downstream in 
order to obtain accurate results. 

In the wing analysis of this paper, an effective wake of 


four root chord lengths was used. The wing program will, 
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however, incorporate any effective wake length desired by 
varying an input parameter (Appendix B). With a wing grid 
Mmmm concurol points, four root chord lengths of effective 
wake requires a wake grid which varies from 400 singularities 
for a rectangular wing, to approximately 1100 singularities 
for a wing with a taper ratio of one-fourth. The increased 
number of wake singularities are required for a tapered wing 
because of the decreasing chordwise dimension of the control 
mexmeoeeowards the wing tip. Thus a greater number of 
Singularities are required, at a spanwise station outboard 
SimeemenrOOL, tO extend four root chord lengths into the eran 
than are required for a rectangular wing. The large number 
Of wake singularities pose no computational problem, since 
the effect of the wake is incorporated prior to the solution 
meee nerne! function matrix equation, as discussed in 


seecumon LV.A.2. 
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B. WING/BODY ANALYSIS 
1. Comparison with Steady Analyses 

A thorough literature search was made in an attempt 
vomurmniad previous theoretical or experimental work, in the 
nonsteady lifting surface field, with which to compare and 
validate the wing/body program results. Unfortunately, no 
such effort which explicitly analyzed body interference 
Seeeeoto On an Oscillating wing's pressure distribution could 
be found. Rodden [38] includes some body interference effects 
jmoears doublet lattice method for analyzing nonplanar con- 
memiieavlons, by extending the lifting surface elements onto 
the body surface near the wing-body intersection. However, 
the actual interference effects are not presented. 

It was therefore necessary to compare the wing/body 
program results with data previously obtained in the steady 
field. The two latest works found were presented at the 
AGARD conference on Aero-dynamic Interference in 1970 by 
Kuchemann [52] and Labrujere [53]. Both analyzed an aspect 
ratio six rectangular wing at a six degree angle of attack 
feeaiommbted On a Cylindrical body of radius approximately 
equal to 20% of the wing semi-span. This configuration was 
also theoretically analyzed by Weber [40] and experimentally 
investigated by Korner [54]. Each of the theoretical analyses 
employed a combination of surface source singularities and 
Tieomeblinc GCistributions to construct a steady kernel func- 
meememet nod solution. - 

Figure 28 compares the wing/body program results, for 


both the wing and wing/body combination, with results 
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obtained from the above investigations. The discrete 
potential element approach spanwise loadings are seen to 
compare almost exactly with the theoretical curves obtained 
by Weber, while the experimental data of Labrujere and 
Korner offer good correlation. The theoretical work of 
Labrujere overestimates the spanwise loading, especially in 
the case of the wing/body combination. The good agreement 
between these theoretical and experimental results and the 
wing/body program results provides a reasonable measure of 
confidence in the discrete potential element analysis of 
Mieomoype Of nonplanar configuration. 

Figure 29 indicates the convergence of the wing/body 
MPieeeram as the number of control points on the quarter cir- 
cumference of the body was increased from (2 wee lO onm Linas 
graph is a greatly expanded representation of the wing root 
area effects. Outboard of the 25% semi-span point, the 
three grid configurations gave essentially the same results 
shown in Figure 28. In each case, body and wing control 
panels were matched as closely as practical in size and 
shape, so that the body grid becomes in essence an extension 
of the wing network on the body surface. Small variations 
between wing and body panel size did not effect the wing 
pressure distribution; however, when wing and body grids 
were obviously mismatched results became inconsistent. 

For all the data presented in this paper, the body grid 


extended from one chord length upstream of the wing leading 


tion point. Extending the body grid somewhat further in 
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either direction produced no appreciable effect on the 
meouevs. However, decreasing the length of the body grid, 
to the point where the wing/wake singularity grid extended 
Deyond the body grid, caused the wing pressure distributions 
Demeceome inconsistent. This effect held for both the 
steady and nonsteady cases. Therefore, it would seem that 
the body surface which effects the wing pressure distribu- 
meen. as Giscussed in Section IV.B.3, is concentrated in 
this effective length between one wing chord upstream and 
downstream of the wing/wake grid. Further investigation of 
this effect, such as by increasing the effective wake and/ 
Greene DOdy grid length with increased number of control 
meno, Was precluded cue to the accuracy limitations of 
Wiemmiaurix inversion routine discussed in Section Vil Ages 

iomtod1l1On, wang body grids within this effective body 
fenmetn were limited in the number of control points which 
could be used for the same reason. It was, therfore, not 
possible to investigate configurations with body radii 
greater than about 20% of the wing semi-span. Future work 
with this method would definitely require an improved com- 
plex matrix inversion procedure to allow greater range in 
the wing/body analysis. 

Mnclcationm Of the effect of body size, in the steady 
case, is presented in Figure 30, where spanwise loadings 
are plotted, with a greatly expanded ordinate, for the wing 
previously considered -in combination with bodies of three 


Gitfrerent radii. lemean be seen that the interference 
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effect is large even for a relatively small body radius of 
five per cent of the wing semi-span. This is reasonable 

Since the body is at zero angle of attack and is providing 
Homelite Carry—-through, or reinforcement, for the wing. This 
effect for bodies of small radius is also shown by Kuchemann's 


results [52]. 
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pemeeexcvensi0n to Nonsteady Interference 


Investigation of wing/body interference for the non- 

Beeady case was conducted for two modes of wing motion: 
bending, as previously discussed; and torsional oscillations 
of the wing cantilevered at the root in an assumed first 
Natural mode. In each case the coefficients were normalized 
with respect to the wing tip angie of attack. It should 
be noted that, of the types of motion normally considered 
in the nonsteady case, these are the only two valid for the 
wing/body configuration considered in this investigation. 

ieee 31 presents the interference effects for the wing 
analyzed in Figure 16, cantilevered from a body with ten per 
cent semi-span radius. The decrease in sectional lift 
amplitude follows the steady case format, being greatest at 
the root and essentially disappearing as the wing tip is 
Pomreached. The change in phase angle, while small, remains 
relatively constant until midspan and then slowly decreases. 
Figure 32 presents spanwise loading in the root area for 
wing/body configurations with body radii varying from five 
to twenty per cent wing semi-span. Again the trends are 
essentially the same as in the steady case, with a small but 
relatively constant phase angle difference. Wing/body inter- 
ference for the torsional vibration case, with body radius 
equal to 20% of the wing semi-span is shown in Figure 33. 
These interference effects in the torsional loading follow 
the same trends as for the nonsteady bending mode, as well 


as the steady case. 
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The results of the wing/body program show that nonsteady 
miveriterence occurs in the same way and via the same mech- 
anism as in the steady case. Interference is greatest in 
the root area, decreasing towards the wing tip. The dif- 
ference between the steady and nonsteady cases comes from 
Bme type of wing loading with which the body interfers. For 
awing at a steady angle of attack, the wing loading is 
greatest in the root area, and, therefore, the interference 
is of relatively large magnitude. In the nonsteady case, 
the wing loading is smaller in the root area (since the wing 
is cantilevered from the body and has no notion at the root) 
and increases towards the tip, as the wing motion amplitude 
increases. Therefore, the interference effects, concentrated 
towards the root area, are of smaller relative magnitude for 


nonsteady motion, as compared to the general steady case. 
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vit. “CONCLUSIONS 


In summary, a discrete potential element approach to 
subsonic numerical lifting surface theory has been developed 
and shown to be practical in predicting the nonsteady load- 
ing on harmonically oscillating wings. This approach was 
then extended to the case of an oscillating wing, canti- 
meer co from a steady cylindrical body, to investigate inter- 
ference effects and show the versatility of the basic method. 

Correlation of the wing program results with those of 
nonsteady lifting surface theory is generally quite good 
Over the full range of subsonic flow. Deviations, where 
they exist, appear to come from the different methods 
employed to handle the wing edge eee re in the pres- 
mee Oistribution, which are assumed a priori in the lifting 
surface theory, but in the discrete potential element 
fPopmoach are implicit in the boundary conditions. Primary 
Geéviation appears to rest in the handling of the wing tip 
pressure slope singularity, since sectional lift and pitching 
moment are generally overestimated in the wing tip area. 
Future work with the discrete potential element method should 
include investigation of a more adequate acknowledgement of 
mise poundary condition in the problem formulation. 

Lmieitimum Of aboultsixty control points is required on 
the semi-span wing surface to achieve convergence of the 


fee Orogram results. In addition, at least six control 
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Mernts are required, in both the chordwise and spanwise 
O@arections, to achieve accurate integration of the pressure 
distribution into sectional and wing coefficients. The 
Mees, which can be analyzed efficiently by this method, are 
effectively limited to medium to small aspect ratios, because 
Meme requirement for large numbers of control points, and 
Ponsequently long computer run times, for high aspect ratio 
fetes.) Lhe wings is further restricted to constant leading 
and trailing edge sweeps, and to a finite tip chord. Future 
development of the program capability should be pointed at 
removing these restrictions. 

The effect of the wake on the wing loading is concen- 
trated in the area just downstream of the wing. In all cases 
investigated, wing loading converged to within one per cent 
dn an effective wake length of four root chords downstream 
of the trailing edge, allowing the wake singularity grid to 
memcerminated at this finite distance from the wing. In 
addition, the wake effect is included in the wing kernel 
mimenton matrix prior to the inversion process, through 
application of the boundary conditions in the wake. In 
this way, an historical drawback to the velocity potential 
formulation, the requirement to integrate the lifting 
surface downwash equation over both wing and wake, is 
iemoved. 

Wing/body program results agree very well with both 
theoretical and experimental results for the steady case. 


4- 


Results obtained for nonsteady wing motion appear consistence 
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and give a good representation of interference effects. 
Convergence of the wing/body program is achieved with a body 
control panel grid of essentially the same format as the 
wing grid. 

The effective body length which causes interference in 
Mme Wing pressure loading extends approximately from one 
chord length upstream of the wing leading edge to one chord 
length downstream of the effective wake termination point. 
Increasing the body length modeled did not effect the 
pressure distribution, while decreasing this length so 
that the body grid did not extend beyond the wing/wake 
eeu larity grid in either the upstream or downstream direc- 
tions caused inconsistent results to occur. This analysis 
feeecomewhnat limited due to numerical restrictions on the 
allowable size of the body grid. 

In the numerical analysis, the number of control points 
on the wing and the body was limited due to loss of compu- 
meaudonal accuracy in the matrix inversion routine of sub- 
routine COMAT. This restricted the scope of the wing/body 
interference investigation both as to the length and radius 
Buebedy which could be considered. Accurate solution of 
the wing downwash equation and inversion of the body influence 
Peemanperent matrix, is limited to coefficient matrices of 
order less than 110. Future wing/body analyses by this 
method would definitely require a more sophisticated inver- 


sion routine capable of handling much larger body grids. 





Within the limitations noted, the wing/body program 
results show that interference effects are significant, and 
follow the same format, in both the steady and nonsteady 
cases. The decrease in wing pressure loading amplitude, 
Gaused by the body's presence, is greatest in the root area, 
decreasing towards the tip. Differences in phase angle, 
while small, exist over the entire wing. Therefore, these 
interference effects are of importance to three-dimensional 
analyses of wing/body configurations. 

Future development of the discrete potential element 
approach should include the analysis of oscillating wings 
with control surfaces. This approach appears to provide a 
merecy means of including the control surface boundary con- 
eictons, without the problems associated with properly 
defining singularity effects in the continuous loading 


formulation. 
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APPENDIX A 


KERNEL FUNCTION DEVELOPMENT 


Hem@eche purposes of nonsteady lifting surface analysis, 
feecetme! function is defined as the normalwash on a surface 
Senco 2 Unit oscillatory acoustic singularity which 
satisfies the linearized convective wave equation (2.9). 

Whe singularities may be elementry radiators of zero order 
representing simple point sources, commonly used in super- 
sonic analyses, or first order radiators, namely dipoles or 
doublets, whose axes are normal to the surface, used in 
SUbsSOnic analyses as in this paper. Development of the 
general form of the doublet kernel function was first 
accomplished by Kussner in 1940 [2]. The Fee lonneat sme 
this appendix follows Kussner's method and specializes the 
Meoults to the forms of the kernel function employed in 
the discrete potential element approach for both the wing 
mecmoody SUuriaces. 

in summary, the solution to the linearized wave equation 
femme sbabionary acoustic singularity is first extended to 
miemease Of a moving singularity in a space fixed coordinate 
Pe eem through application of a Lorentz coordinate transfor- 
mation invariant with respect to the speed of sound. This 
solution is then transferred to the moving wing (or body) 
fixed coordinate system through a Galilean transformation. 


The results are subsequently reduced for the harmonic case 
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imomene iorms of the kernel functions necessary in the analy- 


sis discussed in Section IV. 


ims cabvionary Singularity 
The perturbation analysis of Section II led to the 
linearized convective wave Squavicne (a. switche tor the 


stationary case (U. = 0) has the familiar form 


Vo--—-s —-~= = 0 Ci, 


minee the pressure and velocity potential are linearly 


related by 


p=- po. — Cho) 


Be 
yop - 4 SB = 0 (A.3) 
a ot 


Considering a stationary source at the origin of the coor- 


dinate system, equation (A.3) becomes 


= oa 
1 2 [? 8] ee So 
me Ip se] . — —— = 0 Cres 
re or or aé sto 
where 
ie = xe 32 oe ae Z° 


This equation has the well known solution 


5 = i F (t-r/a,) + 2G (ttr/a,) (A.5) 


i Som 





where F and G are arbitrary functions. From the boundary 
@onditions of the problem, only the solution representing 
Outgoing waves (F) is admitted. 

Representing the oscillating source as a sphere of radius 
Py expanding and contracting with a rate of flow away from 


the surface defined by 


= 2 - 
Set = Wor, u(t) 


The momentum equation (2.2) becomes in linearized form 


3 else ple Be a O oUy, 

ror 2 o dt 
or 

p 
0 co) 6 ols: 
= = ES = ee ae at r= Lh (A.6) 
r unr 
b 
: FF. oF 
if r, is very small, —=z is much larger than ~— at r-r_, and 
b ae OG b 
equation (A.6) can be taken as 
p 
_ 2 Clsiee, _ 

PF in mee ae r 

marendine this solution to a general r 
Po 3 

F = in at S(t-r/a,) 
er 

= P x a 

p> = Unr at See ache), Chey a) 
and 

pos. = (A.7b) 

oO = ie Set -0/ 2.) 1. 1D, 





from equation (A.2). Equations (A.7) represent the pressure 
and velocity potential distributions of a stationary oscil- 


Mavory source of strength 5S. 


moving Singularity 

Consider a source moving with uniform velocity U,, with 
mespect CoO the surrounding fluid, in the direction of the 
negative x axis. If Q(x,y,z,t) is the source distribution 


density, the continuity equation (2.1) becomes 


oP] 
1) 


+ Voll = Q 


Q 


t 


and the linearized wave equation (A.3) becomes 


a O- 
ap — (A. 8) 
2 2 ot 
om ot 


¢ 


mars source distribution can be represented by 


Q(x,y,z,t) = p.S(t) &(xtU,t) sly) 6(z) 


where 6 represents the familiar Dirac delta function. Making 
use of the linear relationship between p and $ in equation 
(A.2), equation (A.8) may then be expressed in terms of the 


velocity potential as 


—~ = S(t) S(x+U,t) d(y) 6(z) (A.9) 


Peierentz transformation, scaled with respect to the compres- 
pipadaty factor 6, can now be used to reduce the right hand 
side of equation (A.9) to that of a stationary source. The 


transformed coordinate system (primed) is defined as follows: 





=s 


x) = (EU. ) 
g° 
ree ¥/76 
(A.10) 
m= 7/8 
U 
t' = a (t+ —= x) 
O 
B a 


meritorming this transformation, equation (A.9) takes the 
Porm 


= 


il il 
—= —s = — S(t')d(x' )6(y')6(z') CA. ay) 
a® gE 1° a 


EFquation (A.9) has therefore been reduced to the form of the 
wave equation for a stationary source in the primed coordi- 


nate system, which from equation (A.7) has the solution 


2 | 
$(r',t') = FA s(tt - r'/a,) (A.12) 


imaecerms Of the space fixed coordinate system, the variables 


of equation (A.12) are 


———- 
ree Ct tt) © + 8° (yo 42°) 


xX ) 


Gale 
i 
oom 
ey 
= 
lg 


Meemeterring the solution (A.12) to the wing fixed coordinate 
system now requires the further Galilean transformation 


x = xtU.t, which produces the final result 


b(r,t) = ee s[t+ (MN x-r) | (A.13) 





i 
3° 


a 
co 
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where 


r= ~/x° + ge (y-+2°) 


All variables are referenced to the wing fixed coordinate 


system. 


See Wing Fixed Periodic Singularities 

If an harmonically oscillating source of unit strength 
(S = eo) is located at the origin of the wing fixed 
eoordinate system, the velocity potential distribution given 


by equation (A.13) is 





= 2 i ; [ a _ i} 
o.(x,y,2,t) = Sie cede {iw ai ie CUR ae) (A.14) 
Considering now only the spatial Varracaon of the velocity 


potential of a unit source located at a general point 


POY 5226) (time dependence eee having been factored out) 





o. (X,Y 52) = Tes exp {3 Fr [™.(x-x,) - R}] ON, 
on 


where 


aos af (x-x,)* + 8° |(y-y,)° + (2-2.)° | 


iimemoerential distribution of a dipole or doublet is 


me lated to that of a source by 


Migemrerm is the direction of the axis of the dipole. There- 


fore from equation (A.15), the velocity potential of a unit 


eimole oriented in the z direction is 


MES 











2. -p° W 
Og0GY 22) 5° (2-2) [1+ 1 2 sfexo(s oft, (xx »- Af 
oo (A.16) 


The downwash at a point on the wing (x,y,0) from such a 


doublet located in the plane of the wing (X59, 20) ese Caer 








by 
J 
d 
w(x,y,0) natal 
WE | iy. yh 
or 
_g¢ 
w= 3 [14 Th 58 exp | R] | Chee) 
KarR aB 
where 


R = ~/(x-x,)° a B(y-y,)° 


This is the standard form of the kernel function in cartesian 
coordinates normally used in the velocity potential formula- 


eromeol! Jifting surface theory. 


4. Kernel Function for Wing/Body Doublets 
To analyze the wing/body problem, a transformation from 
Pomebestan to cylindrical coordinates was made as discussed 


feoection 1V.B and shown in Figure 8. Thus 


Xx = X 
Yas cos 9 
Zar Sin 9 


relate the equations of the previous section to the (r,8,x) 


emer A unit doublet with axis in the z direction located 





aug ea »X) has a potential distribution, from equation 


O 


(A.16), in the new coordinate system of 











2 
—£6 e * ° U) 
ter. o,x) = (Se shaligpmlshe 1a) tealtay 8 [2+ 5] 
d haR? O O a Bo 
exp a [™,,(x-x,) - R} | Oe aalitsy) 


a B° 


where 


7 2 Ce ee 
a= V(x-x,) =F (8) (ie ae = Zio. cos <8-8 >) 


For a doublet located on the wing at (ros 8=0 or T,X); the 


downwash on the wing at (r,o,x) is given by 


Q> 


> 
0 


= 
Tl 
ee Des 
5 


Qo 


v= 0, nae Or 


wim@eem applied to equation (A.18) results in 


i a n)exp {i z i. [M,, (x-x,) = R] | (A.19) 


fo) 


i. = -8° 
e haR? 














where 


ra 
R= (x-x)) + g(r wee + atone 


0 


== Eveiews 


+ for 8 TT 


This is of course the same as equation (A.17) in cartesian 
coordinates except for R which indicates the coordinate 


transformation. 


et 





The normalwash on a cylindrical body of radius Pp» 


Coaxial with the x axis, at a control point, (rp, ,8,x), due 


to the wing doublet is given by 

















a °%a 
ite or a ih 
r=Pp> 8-0 Ore Th 
or 
2 3r 
u, = = z sin of b- e ai faee na 
4atR a8 
o 
+ = rR ee exp {4 st [MC x-x )-R] | (A.20) 
2 B or 2 O 
a,,8 a8 
where 


= O Cy e,.e = 
R= NV (x-x,) + 8B (rptr 4 2rpr. cos @) 


O 
: : 
OR_ B ~ 
ea hl (rp, + rr cos Q) 
- for 0 = 0 
i1o2 0 = TT 


Memdectermine the body singularity kernel function forms, 
maemvelocity potential of a source is expressed in cylindri- 


Pel coordinates. From equation (A.15) 





o.(r,0,x) = — exp {4 Zs. [M,,(x-x,) - | | (A.21) 
ah 


where 


a Oe ee 
aS (x-x,) + B°(r 12m Os Set) 


}— 
ee) 
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mae povential distribution of a doublet oriented in the 


radial direction and located at (r 58 »X) is obtained from 


O 


d. = Oe 
d ar 
O 


which applied to equation (A.21) gives 











> (7, 8,x) = — [ro-r cos (0-0) [2+ 5] 
e xp { 3 2 [t¢x—x,) - R} } (A.22) 
Oo 


With R as above. From this potential formulation, the 
morlowing kernel functions can be obtained for a unit doublet 
with axis in the radial direction, located on the body sur- 
face at (rp, Q 


Ss x5). 


The downwash on the wing at (r,o,x) is given by 











2 OE 
oC) eee 
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Pre Sin 6 
B O 


The normalwash on the body at (rp, ,6,x) is given by 


u 
ie 
or 
u 
it 
where 
R 
aR 
ie 


Balecions (A.19), 


0% 4 
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= rp, (l-cos Se allege, 


(A.20), (A.23), (A.24) represent the 


meme t function formulations used in the discrete potential 


element analysis of the wing and wing/body programs. 





APPENDIX B 


INPUT INSTRUCTIONS FOR COMPUTER PROGRAMS 


meseructions for preparing input data for both the wing 
and wing/body programs are presented in this appendix. The 
Mmeeto location and format for each input quantity is speci- 
fied. Any set of units may be used for geometric dimensions, 
Gisplacements, and acoustic velocity as long as they are 
eonsistent. Any number of problems may be run in sequence, 
with the programs terminating when a new data set is not 


available to be read. 


meeevane Program 
a. First card: FORMAT(3F10.4,2110,F10. 4) 


Colum 1-10 11-20 21-30 LO 50 51-60 


Name XM F A Ce la ALPH 
XM MENGlloly Ja) DuUle eae 

F Circular frequency (rad/sec). 

A Speed of Sound. 

Ga! iidi Cavor slOreseconGsceorceinouy data. 


= 0 for mew data; orogranm reads Second card. 
= 1] for same data as previous problem; no 
Sle(erepael ereecl ive bulieecl- 


JCH Indicator for wing displacement input data. 
= 0 for new data; vrogram reads third and 
subsequent wing displacement cards. 
= ] for same data as previous problem or 
for pirching mMoulLoMmwmmoanine displacement 
cards required. 


ALPH Amplitude of pure pitching motion. 
PEeetOr O15 ChIiie Vereen, Or@eram COmDULeS 
downwash for pitching mode (JCH must equal 1). 
= 0.0 for general wing motion; program com- 
putes downwash from wing displacement input. 
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second card: FORMAT(41I5,5F10.4) 


Column 1-5 6-10 11-15 16-20 21-30 31-40 41-50 51-60 61-70 


Name N M MB NC XC XB BR Gll G2e 
N Number of chordwise control points (maximum ten). 
M Number of spanwise control points (maximum ten). 
MB = 0 for wing program. 

NC Significant wake length downstream from wing 


trailing edge (root chord lengths). 


XC Wing root chord. 

BR = 0 for wing program 

XB Wing semi-span. 

eam Leading edge sweep angle (deg.). 
G22 ga Plaiieeedse. seep wane lem (deg): 


Gil and G22 positive for downstream sweep. 


Third and subsequent wing displacement cards: 
BORAT ( 7F10.4) 


Comumae 1-10 11=20 21-30 37—=40 41=56 51-60 61-70 
Name TG) C2) ECs) ZC eZ (G) 7.C 7) 
Pee) aS) eer at fourth card 
Zl): 2 (le ee Ol 3 Rete iii card 


Displacements (wing motion amplitudes) are read span- 
wise starting from the wing root at the leading edge. 
A new card must be started for each spanline of data 
(example above is for ten spanwise points). 


fee wing/Body Program 


a. 


log 


First card: same as wing program. 
Second card: FORMAT(41I5,3F10,15) 


Column 1-5 6-10 11-15 16-20 21-30 31-40 41-50 51-55 
Name N M MB NC XC XB BR MS 


N Number of wing chordwise control points 
(maximum ten). 


M Number of wing spanwise control points 
(maximum ten). 
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MB 


NC 


XC 
XB 
BR 


MS 


hewa LaMuniber .Ol COMerol Peinuswon body 
quarter-circumference surface (maximum 130) 


Significant wake length downstream from wing 
trailing edge (chord lengths). 


Nal ioe ele\oneral 
Wing semi-span. 
Body radius. 


Number of control points along body quarter 
circumference. 


Third and subsequent wing displacement cards: same 
as wing program. 
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APPENDIX C 


COMPUTER PROGRAM FLOW DIAGRAMS 
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FIGURE 36 
_DINCO FLOW DIAGKAM 
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FIGURE 37 


SUBROUTINE COMAT FLOW DIAGRAM 
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FIGURE 38 


SUBROUTINES SECLM AND WINGLM FLOW DIAGKAM 
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FIGURE 39° 
\AIYNG/BODY PROGRAM - MAIN FLOW DIAGRAM 
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FIGURE 40 


WING/BODY PROGRAM - DINCO FLOW DIAGRAM 
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SUBROUTINE UTHETA FLOW DIAGRAM 
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FIGURE 472 
SUBROUTINE URADL FLOW DIAGKAM 
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FIGURE 43 
SUBROUTINE URKADZ FLOW DIAGKAM 
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